In [69]:
import pandas as pd
import numpy as np
from scipy.integrate import trapz

import importlib
import useful
importlib.reload(useful)

import matplotlib
import matplotlib.pyplot as plt
#parameter to use Latex in matplotlib
matplotlib.rcParams.update(matplotlib.rcParamsDefault)
#a plot label should be easy to read
plt.rcParams.update({'font.size': 14})
import matplotlib.image as mpimg

Cospectrum Budget Model Analysis of FLOSS data

Lorenzo Lazzarino supervised by Luca Mortarini

  • The focus here is to study spectra and cospectra of turbulent wind motion (u,v,w) over a snow-covered surface at several height levels.
  • FLOSS project: studying the surface meteorology of snow-covered rangeland in the North Park region of Colorado, near Walden (https://data.eol.ucar.edu/project/FLOSS-II)
  • the final aim of this short project is to apply the Cospectrum Budget Model to predict the $F_{uw}$ cospectrum from the (papers in pdf are in this folder), namely:
In [2]:
# Load and display the image
img = mpimg.imread("cospectrumModel.png.png")
plt.imshow(img)
plt.axis("off")  # Hide axes
plt.show()
  • notably, this complete model encompasses also the fluxes of u and temperature. In a simplified version, we will use the model without the $F_{u,\theta}$ term
  • the model here is written in terms of the wave number $k_x$, defined as $2 \pi * \text{frequency}/\bar{u_x}$, where $\bar{u_x}$ is the mean wind velocity along the x direction at the particular height considered
Briefly on turbulence:
  • By analyzing the velocity of wind (let's say u, along the gradient of pressure and parallel to the floor) in the frequency (and wavenumber) domain, we can get a clear picture of scales.
  • Low-frequency scales can be regarded as submeso motion, waves, and turbulence in the form of big, slow vortices that eventually break up into smaller vortices at higher frequencies.
  • This implies a transfer of energy from low-frequency scales to high-frequency scales, where no energy is produced or dissipated.
  • We expect to observe a peak in the power spectrum of velocity fields where energy is produced, followed by a cascade. towards higher frequencies. This behaviour is better identified in w spectrum since it is less affected by the large scale.
  • The cascade's behavior, as a consequence of Kolmogorov's theory of turbulence, shows a hallmark (-5/3) power law for the components of the wind velocity in the spectrum and (-7/3) in the cospectrum
  • As we analyze turbulence at different height levels, we expect it to be less energetic at higher heights since friction with the ground becomes less relevant.
Algorithms and functions are in useful.py!!!.
In [3]:
import matplotlib.image as mpimg

# Load and display the image
img = mpimg.imread("flowchart.jpg")
# Create a larger figure with higher DPI
plt.figure(figsize=(12, 8), dpi=200)  # Increase DPI and figure size
plt.imshow(img)
plt.axis("off")  # Hide axes
plt.show()

Prepare the data

Read one day of data

In [4]:
df = useful.read_day("apr01.dat")

    u_1   v_1   w_1     T_1   u_2   v_2   w_2     T_2   u_5   v_5  ...  w_15  \
0  5.14  1.39  0.10  277.72  5.12  1.55 -0.06  278.06  6.24  2.02  ... -0.31   
1  5.02  1.30  0.13  277.68  5.20  1.96  0.19  278.12  6.41  1.93  ... -0.52   
2  5.00  1.30 -0.01  277.72  4.90  1.67 -0.08  278.04  6.53  1.76  ... -0.60   
3  4.87  1.26 -0.10  277.72  4.72  1.88 -0.05  278.07  6.50  1.75  ... -0.60   
4  5.01  1.09 -0.16  277.71  4.64  1.82  0.02  278.13  6.67  1.49  ... -0.63   

     T_15  u_20  v_20  w_20    T_20  u_30  v_30  w_30    T_30  
0  280.23  9.47  3.97 -0.20  280.24  9.68  4.47  0.01  280.72  
1  280.11  9.22  3.90 -0.08  280.27  9.65  4.42  0.00  280.71  
2  280.17  9.21  4.00 -0.09  280.27  9.48  4.25 -0.08  280.57  
3  280.17  9.22  3.91 -0.13  280.32  9.66  4.40 -0.02  280.67  
4  280.15  9.26  3.99 -0.10  280.27  9.55  4.30 -0.05  280.74  

[5 rows x 28 columns]

NB: in dataframe df i don't have time (but can use data number) - so I modified the rotation algorithm not to need it

per hour it's 7 tables of u,v,w,T to which apply rotations (it's the seven heights) [an hour is 2.16e5 rows]

Sketch of the analysis

In [5]:
# separate in 10 hours


# separate in 7 heights


# on one height: 
#   - plot u,v,w (just to start)
#   - rotate hourly
#   - test the rotation
#   - plot u,v,w,T, \sigma^2_u, \sigma^2_w, \sigma^2_uw, \sigma^2_wT #kinetcik energy components and fluxes (u,w)

# second part:
#   - spectra and cospectra
#   - sigma squared from the integral
#   - smoothing spectra
#   - nice plotting
#   - apply model

split data in matrices

In [6]:
df_matrix = useful.split_dataframe(df, 10, 7)
In [7]:
df_matrix[0][0]
Out[7]:
u_1 v_1 w_1 T_1
0 5.14 1.39 0.10 277.72
1 5.02 1.30 0.13 277.68
2 5.00 1.30 -0.01 277.72
3 4.87 1.26 -0.10 277.72
4 5.01 1.09 -0.16 277.71
... ... ... ... ...
215995 3.52 2.69 0.55 276.44
215996 3.66 2.73 0.37 276.34
215997 3.55 2.82 0.41 276.46
215998 4.13 2.53 0.07 276.78
215999 3.82 2.81 0.23 276.55

216000 rows × 4 columns

split data in dictionaries (will be using this)

  • one day of aggregated data will be split in hours and heigths
In [8]:
row_keys = ["h20","h21","h22","h23","h0","h1","h2","h3","h4","h5"] #hour of measurement
col_keys = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"] #livelli FLOSS
df_matrix_dict = useful.split_dataframe_dict(df, 10, 7, row_keys, col_keys)
In [9]:
df_matrix_dict["h5_z30"] 

# in (i,j): f"{row_keys[i]}_{col_keys[j]}" to iterate
# or for key, chunk in df_dict.items():
#    print(f"\nDataFrame {key}:")
#    print(chunk)
Out[9]:
u_30 v_30 w_30 T_30
1944000 -0.93 6.55 -0.03 276.57
1944001 -0.96 6.57 -0.01 276.56
1944002 -0.96 6.57 -0.03 276.56
1944003 -0.94 6.59 -0.03 276.58
1944004 -1.02 6.57 -0.03 276.59
... ... ... ... ...
2159995 -3.85 5.31 0.31 275.02
2159996 -4.01 5.31 0.31 275.04
2159997 -3.85 5.40 0.30 275.03
2159998 -3.77 5.48 0.29 275.00
2159999 -3.86 5.29 0.23 274.97

216000 rows × 4 columns

let's take a look at the data

In [10]:
useful.plotHour(df_matrix_dict, hour= "h5", height = "z30")

Rotation

we need to align (x,y,z) and thus (u,v,w) to perpendicular to the floor and along the local hourly pressure gradient and we do so by applying three rotation with the requests:

  • $\langle v \rangle = 0$
  • $\langle w \rangle = 0$
  • $\langle v'w' \rangle = 0$
In [11]:
df_matrix_dict["h5_z30"] = useful.rotateHour(df_matrix_dict, hour= "h5", height = "z30")
In [12]:
df_matrix_dict["h5_z30"]
Out[12]:
u_30 v_30 w_30 T_30
1944000 6.562963 -0.833803 -0.024334 276.57
1944001 6.590320 -0.808850 -0.006198 276.56
1944002 6.590165 -0.810208 -0.026151 276.56
1944003 6.604183 -0.834725 -0.024591 276.58
1944004 6.605987 -0.752474 -0.030203 276.59
... ... ... ... ...
2159995 6.139542 2.325891 0.104702 275.02
2159996 6.181734 2.479849 0.093895 275.04
2159997 6.226276 2.301487 0.095665 275.03
2159998 6.282268 2.202740 0.091927 275.00
2159999 6.122267 2.335353 0.024004 274.97

216000 rows × 4 columns

Analysis

u,v,w,T rotated time series

In [13]:
useful.plotHour(df_matrix_dict, hour= "h5", height = "z30")
In [14]:
useful.plotHour(df_matrix_dict, hour= "h20", height = "z30")

Vertical profile (and gradient) of umean

Let's take a look at the horizontal mean u velocity at each height for the hours h5 and h20. The gradient of this vertical profile will be necessary for the cospectrum budget model.

In [15]:
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
hour="h5"
umeans = []
for height in heights:
    hour_loc = hour + "_" + height
    df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
    u = df_matrix_dict[hour_loc].iloc[:,0].values
    umeans.append(np.mean(u))
In [16]:
heightsmeters = np.array([1,2,5,10,15,20,30])
# Creazione del grafico
plt.figure(figsize=(4,3))
plt.scatter(umeans,heightsmeters, marker='o', color='coral')

# Miglioramenti estetici
plt.title(r"$\bar{u}(z)$")
plt.ylabel("Height $[m]$")
plt.xlabel(r"$\bar{u} \, [m/s]$")
plt.grid(True, which='both', linestyle='--', linewidth=0.5)

# Mostra il grafico
plt.show()
In [17]:
from scipy.interpolate import CubicSpline

# Create the cubic spline
cs = CubicSpline(heightsmeters,umeans)

# Generate a smooth set of x values for plotting the spline fit
x_fit = np.linspace(heightsmeters.min(), heightsmeters.max(), 200)
y_fit = cs(x_fit)

plt.figure(figsize=(4,3))
plt.scatter(umeans,heightsmeters, marker='o', color='coral')

plt.title(r"$\bar{u}(z)$")
plt.ylabel("Height $[m]$")
plt.xlabel(r"$\bar{u} \, [m/s]$")
plt.grid(True, which='both', linestyle='--', linewidth=0.5)

plt.plot(y_fit, x_fit, label="Cubic Spline Fit", color='black', linewidth=2)


plt.legend()

plt.show()
  • clearly we need some smoothing in the fit ...

https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CubicSpline.html

In [18]:
from scipy.interpolate import UnivariateSpline

# Create a smoothing spline with a chosen smoothing factor (s)
smoothing_factor = 0.5  # Adjust this for more or less smoothing
spline = UnivariateSpline(heightsmeters,umeans, s=smoothing_factor)

y_fit = spline(x_fit)  # Smoothed spline

plt.figure(figsize=(4,3))
plt.scatter(umeans,heightsmeters, marker='o', color='coral')

plt.title(r"$\bar{u}(z)$")
plt.ylabel("Height $[m]$")
plt.xlabel(r"$\bar{u} \, [m/s]$")
plt.grid(True, which='both', linestyle='--', linewidth=0.5)

plt.plot(y_fit, x_fit, label="Smoothed Spline Fit", color='black', linewidth=2)


plt.legend()

plt.show()
In [19]:
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
hour="h20"
umeans = []
for height in heights:
    hour_loc = hour + "_" + height
    df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
    u = df_matrix_dict[hour_loc].iloc[:,0].values
    umeans.append(np.mean(u))
    
from scipy.interpolate import UnivariateSpline

# Create a smoothing spline with a chosen smoothing factor (s)
smoothing_factor = 0.5  # Adjust this for more or less smoothing
spline = UnivariateSpline(heightsmeters,umeans, s=smoothing_factor)

y_fit = spline(x_fit)  # Smoothed spline

plt.figure(figsize=(4,3))
plt.scatter(umeans,heightsmeters, marker='o', color='coral')

plt.title(r"$\bar{u}(z)$")
plt.ylabel("Height $[m]$")
plt.xlabel(r"$\bar{u} \, [m/s]$")
plt.grid(True, which='both', linestyle='--', linewidth=0.5)

plt.plot(y_fit, x_fit, label="Smoothed Spline Fit", color='black', linewidth=2)


plt.legend()

plt.show()

Vertical profile of mean Potential Temperature

it is interesting also to take a look at the potential temperature profile $\theta_v = 0.0098*\bar{T}*\text{height}$

In [20]:
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
hour="h5"
Tmeans = []
for height in heights:
    hour_loc = hour + "_" + height
    df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
    T = df_matrix_dict[hour_loc].iloc[:,3].values
    Tmeans.append(np.mean(T))
    
Tmeans = np.array(Tmeans)
heightsmeters = np.array([1,2,5,10,15,20,30])
# Creazione del grafico
plt.figure(figsize=(4,3))
plt.scatter(Tmeans*0.0098*heightsmeters,heightsmeters, marker='o', color='rebeccapurple')

# Miglioramenti estetici
plt.title(r"$\theta_v(z)$, surface is snow")
plt.ylabel("Height $[m]$")
plt.xlabel(r"$\theta_v(z) \, [K\cdot m]$")
#plt.xlim(270,280)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)

# Mostra il grafico
plt.show()
In [21]:
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
hour="h20"
Tmeans = []
for height in heights:
    hour_loc = hour + "_" + height
    df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
    T = df_matrix_dict[hour_loc].iloc[:,3].values
    Tmeans.append(np.mean(T))

Tmeans = np.array(Tmeans)
heightsmeters = np.array([1,2,5,10,15,20,30])
# Creazione del grafico
plt.figure(figsize=(4,3))
plt.scatter(Tmeans*0.0098*heightsmeters,heightsmeters, marker='o', color='rebeccapurple')

# Miglioramenti estetici
plt.title(r"$\theta_v(z)$, surface is snow")
plt.ylabel("Height $[m]$")
plt.xlabel(r"$\theta_v(z) \, [K\cdot m]$")
#plt.xlim(270,280)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)

# Mostra il grafico
plt.show()
  • temperature after nighttime (h5) is lower
  • also right after nighttime the atmosphere is more stable (usualy weakly and very stable stratification alternate at nighttime)

Covariance

Before calculating spectra and cospectra, one can evaluate sigmas and covariances, which represents turbulent energies and turbulent fluxes, which are terms of the turbulent kinetic energy equation.

  • for a specific height:
In [22]:
covmat = np.cov(df_matrix_dict["h5_z30"].values,rowvar=False)
print(covmat)

[[ 6.51908258e-01  2.36527488e-01 -5.03899950e-02  3.78975550e-01]
 [ 2.36527488e-01  5.04895780e-01 -5.33731910e-18 -4.96202925e-02]
 [-5.03899950e-02 -5.33731910e-18  5.46453484e-02 -5.84562002e-02]
 [ 3.78975550e-01 -4.96202925e-02 -5.84562002e-02  6.85336962e-01]]
In [23]:
from IPython.display import display, Markdown

a = 1.23
b = 4.56
c = 7.89
d = 0.12

table_latex = f"""
\\
\\begin{{array}}{{|c|c|}}
\\hline
\\textbf{{energia/flussi}} & \\textbf{{valori}} \\\\
\\hline
\\text{{$\sigma^2$$_u$$_u$}} & {covmat[0,0]:.4f} \\\\
\\hline
\\text{{$\sigma^2$$_w$$_w$}} & {covmat[2,2]:.4f} \\\\
\\hline
\\text{{$\sigma^2$$_u$$_w$}} & {covmat[0,2]:.4f} \\\\
\\hline
\\text{{$\sigma^2$$_w$$_T$}} & {covmat[2,3]:.4f} \\\\
\\hline
\\end{{array}}
\
"""

display(Markdown(table_latex))

\ \begin{array}{|c|c|} \hline \textbf{energia/flussi} & \textbf{valori} \\ \hline \text{$\sigma^2$$_u$$_u$} & 0.6519 \\ \hline \text{$\sigma^2$$_w$$_w$} & 0.0546 \\ \hline \text{$\sigma^2$$_u$$_w$} & -0.0504 \\ \hline \text{$\sigma^2$$_w$$_T$} & -0.0585 \\ \hline \end{array}

  • profiles:
In [24]:
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
hour="h20"
suu = []
sww = []
suw = []
swT = []
for height in heights:
    hour_loc = hour + "_" + height
    df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
    covmat = np.cov(df_matrix_dict[hour_loc].values,rowvar=False)
    suu.append(covmat[0,0])
    sww.append(covmat[2,2])
    suw.append(covmat[0,2])
    swT.append(covmat[2,3])
    
heightsmeters = np.array([1,2,5,10,15,20,30])
# Set up a 2x2 grid of subplots
fig, axs = plt.subplots(2, 2, figsize=(7,6))

# Plot each quantity on a separate subplot
# First subplot for `suu`
axs[0, 0].plot(suu, heightsmeters, marker='o', color='blue')
axs[0, 0].set_xlabel(r"$\sigma^2_{uu} [m^2/s^2]$")
axs[0, 0].set_ylabel("Height $[m]$")
axs[0, 0].grid(True, which='both', linestyle='--', linewidth=0.5)

# Second subplot for `sww`
axs[0, 1].plot(sww, heightsmeters, marker='o', color='green')
axs[0, 1].set_xlabel(r"$\sigma^2_{ww} [m^2/s^2]$")
axs[0, 1].grid(True, which='both', linestyle='--', linewidth=0.5)

# Third subplot for `suw`
axs[1, 0].plot(suw, heightsmeters, marker='o', color='red')
axs[1, 0].set_xlabel(r"$\overline{u'w'} [m^2/s^2]$")
axs[1, 0].set_ylabel("Height $[m]$")
axs[1, 0].grid(True, which='both', linestyle='--', linewidth=0.5)

# Fourth subplot for `swT`
axs[1, 1].plot(swT, heightsmeters, marker='o', color='purple')
axs[1, 1].set_xlabel(r"$\overline{w'T'} [m^2/s^2]$")
axs[1, 1].grid(True, which='both', linestyle='--', linewidth=0.5)

# Adjust layout and show plot
plt.tight_layout()
plt.suptitle("Covariance components vs Height", y=1.02, fontsize=14)
plt.show()
In [25]:
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
hour="h5"
suu = []
sww = []
suw = []
swT = []
for height in heights:
    hour_loc = hour + "_" + height
    df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
    covmat = np.cov(df_matrix_dict[hour_loc].values,rowvar=False)
    suu.append(covmat[0,0])
    sww.append(covmat[2,2])
    suw.append(covmat[0,2])
    swT.append(covmat[2,3])
    
heightsmeters = np.array([1,2,5,10,15,20,30])
# Set up a 2x2 grid of subplots
fig, axs = plt.subplots(2, 2, figsize=(7,6))

# Plot each quantity on a separate subplot
# First subplot for `suu`
axs[0, 0].plot(suu, heightsmeters, marker='o', color='blue')
axs[0, 0].set_xlabel(r"$\sigma^2_{uu} [m^2/s^2]$")
axs[0, 0].set_ylabel("Height $[m]$")
axs[0, 0].grid(True, which='both', linestyle='--', linewidth=0.5)

# Second subplot for `sww`
axs[0, 1].plot(sww, heightsmeters, marker='o', color='green')
axs[0, 1].set_xlabel(r"$\sigma^2_{ww} [m^2/s^2]$")
axs[0, 1].grid(True, which='both', linestyle='--', linewidth=0.5)

# Third subplot for `suw`
axs[1, 0].plot(suw, heightsmeters, marker='o', color='red')
axs[1, 0].set_xlabel(r"$\overline{u'w'} [m^2/s^2]$")
axs[1, 0].set_ylabel("Height $[m]$")
axs[1, 0].grid(True, which='both', linestyle='--', linewidth=0.5)

# Fourth subplot for `swT`
axs[1, 1].plot(swT, heightsmeters, marker='o', color='purple')
axs[1, 1].set_xlabel(r"$\overline{w'T'} [m^2/s^2]$")
axs[1, 1].grid(True, which='both', linestyle='--', linewidth=0.5)

# Adjust layout and show plot
plt.tight_layout()
plt.suptitle("Covariance components vs Height", y=1.02, fontsize=14)
plt.show()

Spectra \& cospectra

I need an algorithm that

  1. performs dft on u,w,T

Then go on spectra/cospectra

  1. evaluates E{uu}, E{ww}, F{uw}, F{wT} of the dfts
  2. plots the spectral functions E
  3. fits the -5/3 trend on high frequencies (and maybe trends on T)
  4. integrates them to get the sigmas back (which percentage, shall find only half of sigma since you ditch left part, the negative one) -> a strategy can be implementing trapezioids integration

Questions:

  • -5/3 "fit"?
  • normalization of mean and variance of signal before dft-ing?
  • which algorithm of dft
  • premultiplied? scale to display?
In [26]:
useful.plotPowerSpectra(df_matrix_dict, hour= "h5", height = "z30", normalize = True,)
In [27]:
useful.plotPowerSpectraPremultiplied_Scatter(df_matrix_dict, hour= "h5", height = "z30", normalize = True)
In [28]:
useful.plotPowerSpectraPremultiplied(df_matrix_dict, hour= "h5", height = "z30", normalize = True)

On high frequencies one can see white noise! it is equidistributed along all frequencies, therefore it is a line with $m=1$ in the premultiplied spectrum

dft algorithm used:

https://numpy.org/devdocs/reference/generated/numpy.fft.fft.html

https://numpy.org/devdocs/reference/generated/numpy.fft.fftfreq.html

the premultiplied spectrum is an energy. its "peak":

  • is in the region of turbulence formation in w
  • is in low frequencies in T

kolmogorov idea: one puts itself in the "reference frame of turbulence"

the integral under the curve

  • integral of E with trapz recoups $\sigma^2$. For example let's take the case of $E_{uu}$, the integral is $0.651905149$ while $\sigma^2$ from the time series is $0.651905240$

The next 4 cospectra plots report only the negative part of the cospectra

since it is the negative part we are interested on (usually one normalizes by the integral/covariance, which is negative)

In [31]:
useful.plotAllCospectra(df_matrix_dict, hour= "h5", height = "z30")
In [37]:
useful.plotAllPremultipliedCospectra(df_matrix_dict, hour= "h5", height = "z30")

Remove noise at high frequency

  • first we do this with a moving average on freq plot. We tried both overlapping and non overlapping windows. The second way was best to handle medium frequencies.
  • the other option is splitting the time series, taking multiples spectra and average them (not done here)
In [38]:
useful.plotPowerSpectra(df_matrix_dict, hour= "h5", height = "z30", normalize = True, movingAver=True, window_size=50, high_freq = 2*10**-2)
In [39]:
useful.plotPowerSpectraPremultiplied(df_matrix_dict, hour= "h5", height = "z30", normalize = True, movingAver=True, window_size=50, high_freq = 2* 10**-2)

Plot smoothed over original spectrum

In [40]:
useful.plotAllCospectra(df_matrix_dict, hour= "h5", height = "z30",movingAver=True, window_size=50, high_freq = 2*10**-2)
In [41]:
useful.plotAllPremultipliedCospectra(df_matrix_dict, hour= "h5", height = "z30",movingAver=True, window_size=50, high_freq = 2*10**-2)

uw is energy extraction, we expect it to be a negative quantity -> since I want to focus on the turbulence part, I will divide the cospectrum by its integral, which is negative for uw

To apply the model we need:

  • uw, uT, ww -> the whole spectra

  • (uu) -> epsilon from a fit with the hallmark (-5/3) power law exponent of turbulence in the velocity fields.

In order to normalize spectra one can use covariance or the proper integral -> I did the latter

now we report both negative (before dividing by integral) and positive part

In [43]:
useful.plotSomeCospectraNegAndPos(df_matrix_dict, hour= "h5", height = "z30",movingAver=True, window_size=50, high_freq = 2*10**-2)
In [45]:
useful.plotSomePremCospectraNegAndPos(df_matrix_dict, hour= "h5", height = "z30",movingAver=True, window_size=50, high_freq = 2*10**-2)

Dissipation of the Turbulent Kinetic Energy

  • fit area: intuitively and as an educated guess, then many moving windows with the same number of points so that each fit has the same statistical features
  • 500 points, overlapping moving windows, since the amount of data is rich

from the f*E_uu plot I can see that a -large- window may be chosen between the frequencies $10^{-1}$ e $10^1$. How many data points are there?

In [46]:
sampling_rate = 60     # 60 Hz
T = 3600               # 1hour duration
N = int(T * sampling_rate)  # Number of samples
frequencies = np.fft.fftfreq(N, 1/sampling_rate) #freq 0 is the mean!
In [47]:
np.sum((frequencies>=0.1) & (frequencies<=10))
Out[47]:
35641
In [48]:
#create a mask
target_freqs = (frequencies>= 0.1) & (frequencies<= 10)
In [49]:
frequencies[target_freqs].size
Out[49]:
35641
In [60]:
#calc E_uu
freqs, powerSpectrumSmooth = useful.CalcPowerSpectrum(df_matrix_dict, hour= "h5", height = "z30", normalize = True,
                                               quantity='u', movingAve=True, windowsize=50, highfreq = 2*10**-2)
fig,ax = plt.subplots(1,1,figsize=(6,3))
ax.loglog(freqs,powerSpectrumSmooth, color="rebeccapurple", marker='o', linestyle='none', markersize=2, ) 
slope = -5/3
y_intercept = 1  # you can change this if you want a different intercept
x_line = np.linspace(0.1, 30, 100)  # line x values
y_line =  x_line**slope* y_intercept
ax.plot(x_line, y_line, linestyle='--', color='red', label=r'$y \sim x^{-5/3}$ ')
ax.grid(True)
ax.legend()
ax.set_ylim(1e-7, 1e4) 
ax.set_xlabel("frequencies [Hz]")
ax.set_ylabel(r"$E_{uu} [m^2/s^2]$")
Out[60]:
Text(0, 0.5, '$E_{uu} [m^2/s^2]$')
In [51]:
#now we fit on moving logarithmically spaced windows
from scipy.stats import linregress

# Reset lists if re-running code
slopes, intercepts, r_values = [], [], []

window_size= 500
step = 1

target_freqs = (freqs>= 0.1) & (freqs<= 10)
freqs_interval=freqs[target_freqs]
spectrum_interval=powerSpectrumSmooth[target_freqs]

for i in range(0, freqs_interval.size - window_size + 1, step):
    window_spec = spectrum_interval[i:i + window_size]
    window_freqs = freqs_interval[i:i + window_size]
    
    # Perform linear regression on each window
    result = linregress(np.log10(window_freqs), np.log10(window_spec))
    
    # Store slope, intercept, and R-squared
    slopes.append(result.slope)
    intercepts.append(result.intercept)
    r_values.append(result.rvalue ** 2)  # R-squared
    
intercepts=np.array(intercepts)
slopes = np.array(slopes)
In [52]:
# Display results
# Create a figure with two subplots for the histograms
fig, axs = plt.subplots(1, 2, figsize=(10, 4))

# Plot histogram for slopes
axs[0].hist(slopes, bins=50, color='blue', edgecolor='black')
axs[0].set_title('Histogram of Slopes')
axs[0].set_xlabel('Slope')
axs[0].set_ylabel('Frequency')


# Plot histogram for intercepts
axs[1].hist(intercepts, bins=50, color='green', edgecolor='black')
axs[1].set_title('Histogram of Intercepts')
axs[1].set_xlabel('Intercept')
axs[1].set_ylabel('Frequency')

# Display the histograms
plt.tight_layout()
plt.show()

I choose the fit with "the best -5/3" (square the values then argmin),

then i build up the line to plot it,

and extract the epsilon

In [53]:
# Display results
# Create a figure with two subplots for the histograms
fig, axs = plt.subplots(1, 2, figsize=(10, 4))
Pslopes = (slopes + 5/3)**2
# Plot histogram for slopes
axs[0].hist(Pslopes, bins=50, color='blue', edgecolor='black')
axs[0].set_title('Histogram of Slopes')
axs[0].set_xlabel('Slope')
axs[0].set_ylabel('Frequency')

# Plot histogram for intercepts
axs[1].hist(intercepts, bins=50, color='green', edgecolor='black')
axs[1].set_title('Histogram of Intercepts')
axs[1].set_xlabel('Intercept')
axs[1].set_ylabel('Frequency')

# Display the histograms
plt.tight_layout()
plt.show()
In [54]:
print(slopes[np.argmin(Pslopes)], intercepts[np.argmin(Pslopes)])

-1.630259334553686 -2.1771024175707634
In [59]:
#calc E_uu
freqs, powerSpectrumSmooth = useful.CalcPowerSpectrum(df_matrix_dict, hour= "h5", height = "z30", normalize = True,
                                               quantity='u', movingAve=True, windowsize=100, highfreq = 10**-2)
fig,ax = plt.subplots(1,1,figsize=(6,3))
ax.loglog(freqs,powerSpectrumSmooth, color="rebeccapurple", marker='o', linestyle='none', markersize=2, ) 
slope = slopes[np.argmin(Pslopes)]
y_intercept = 10**intercepts[np.argmin(Pslopes)]  # you can change this if you want a different intercept
x_line = np.linspace(0.1, 30, 100)  # line x values
y_line =  x_line**slope* y_intercept
ax.plot(x_line, y_line, linestyle='--', color='red', label=r'$y \sim x^{-5/3}$ ')
ax.grid(True)
ax.legend()
ax.set_ylim(1e-7, 1e4) 
ax.set_xlabel("frequencies [Hz]")
ax.set_ylabel(r"$E_{uu} [m^2/s^2]$")
Out[59]:
Text(0, 0.5, '$E_{uu} [m^2/s^2]$')

the model for epsilon is:

$ E_\varepsilon(f)=A \,\bar{u}^{2/3}\varepsilon^{2/3}f^{-5/3} $

with $A=0.15$ for $u$ and $0.2$ for $v,w$

"10**intercepts[np.argmin(Pslopes)]" is $f^{-5/3}$ prefactor

In [64]:
u = df_matrix_dict["h5_z30"].iloc[:,0].values
np.mean(u)
Out[64]:
5.598593989241382
In [65]:
epsilon = (  10**intercepts[np.argmin(Pslopes)] /0.15  )**(3/2) / np.mean(u)
epsilon
Out[65]:
0.0016677483187532535

now plot the spectrum from the epsilon model (I put exactly -5/3)

In [66]:
#calc E_uu
freqs, powerSpectrumSmooth = useful.CalcPowerSpectrum(df_matrix_dict, hour= "h5", height = "z30", normalize = True,
                                               quantity='u', movingAve=True, windowsize=100, highfreq = 10**-2)
fig,ax = plt.subplots(1,1,figsize=(6,3))
ax.loglog(freqs,powerSpectrumSmooth, color="rebeccapurple", marker='o', linestyle='none', markersize=2, ) 

slope = -5/3
x_line = np.linspace(0.1, 30, 100)  # line x values
y_line =  x_line**slope* 0.15 * np.mean(u)**(2/3) * epsilon**(2/3)

ax.plot(x_line, y_line, linestyle='--', color='orange', label=r'$E_{uu} = A \,\bar{u}^{2/3} \epsilon^{2/3} f^{-5/3}$ ')
ax.grid(True)
ax.legend()
ax.set_ylim(1e-7, 1e4) 
ax.set_xlabel("frequencies [Hz]")
ax.set_ylabel(r"$E_{uu} [m^2/s^2]$")
Out[66]:
Text(0, 0.5, '$E_{uu} [m^2/s^2]$')

Let's get some epsilons from u,v,w and all the heights

In [70]:
eps = useful.getEpsilon(df_matrix_dict, hour= "h5", height = "z30", plot=True)

[0.0016677483187532535, 0.0015087758866000156, 0.0009482942703985655]
In [71]:
hours = ["h20","h21","h22","h23","h0","h1","h2","h3","h4","h5"] #hour of measurement
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"] #livelli FLOSS
In [72]:
hour="h5"
epsilons = []
for height in heights:
    hour_loc = hour + "_" + height
    df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
    eps = useful.getEpsilon(df_matrix_dict, hour, height , plot=False)
    epsilons.append(eps)

[0.007399051060660499, 0.007889833127608787, 0.005093578944788387]
[0.006645618167514078, 0.007100642534677185, 0.004423474863920642]
[0.005340602719809727, 0.006076480216627502, 0.003583618004965432]
[0.002454760974526618, 0.002648793160847895, 0.0016084229094623012]
[0.0019572927047798077, 0.0021133558271084684, 0.0013544123909550651]
[0.0017894594821892458, 0.0017828179986896686, 0.001207545039661742]
[0.0016677483187532535, 0.0015087758866000156, 0.0009482942703985655]
In [73]:
hour="h20"
epsilons = []
for height in heights:
    hour_loc = hour + "_" + height
    df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
    eps = useful.getEpsilon(df_matrix_dict, hour, height , plot=False)
    epsilons.append(eps)

[0.06577826542237629, 0.05073024224842888, 0.04573471488874804]
[0.04129473340960273, 0.037393911210352945, 0.03145324377875913]
[0.029444981954559832, 0.028835444383382102, 0.0242597551771929]
[0.018261855843780227, 0.018302642359700583, 0.014910673064156757]
[0.012478869368254451, 0.012840306696191888, 0.01114578570488523]
[0.008523239621536182, 0.00843650217244339, 0.0067336521957292]
[0.0023345902251810675, 0.0019842642782632487, 0.001726669732759961]
In [74]:
# Dati
data = [
    [0.06577826542237629, 0.05073024224842888, 0.045734714888747974],
[0.04129473340960273, 0.037393911210352945, 0.03145324377875913],
[0.029444981954559846, 0.028835444383382078, 0.024259755177192876],
[0.018261855843780227, 0.01830264235970057, 0.014910673064156766],
[0.012478869368254451, 0.012840306696191888, 0.01114578570488523],
[0.008523239621536182, 0.00843650217244339, 0.0067336521957292],
[0.0023345902251810675, 0.0019842642782632513, 0.001726669732759961],
]

# Etichette delle righe e intestazioni delle colonne
row_labels = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
column_labels = ["u", "v", "w"]

# Genera il codice LaTeX per la tabella
latex_code = "\\begin{array}{lccc}\n"
latex_code += " & " + " & ".join(column_labels) + " \\\\\n"
latex_code += "\\hline\n"
for row_label, row_data in zip(row_labels, data):
    row = f"{row_label} & " + " & ".join(f"{x:.6f}" for x in row_data) + " \\\\\n"
    latex_code += row
latex_code += "\\end{array}"

print(latex_code)

\begin{array}{lccc}
 & u & v & w \\
\hline
z1 & 0.065778 & 0.050730 & 0.045735 \\
z2 & 0.041295 & 0.037394 & 0.031453 \\
z5 & 0.029445 & 0.028835 & 0.024260 \\
z10 & 0.018262 & 0.018303 & 0.014911 \\
z15 & 0.012479 & 0.012840 & 0.011146 \\
z20 & 0.008523 & 0.008437 & 0.006734 \\
z30 & 0.002335 & 0.001984 & 0.001727 \\
\end{array}

we now report the epsilon for h5 and h20 in tables

In [75]:
from IPython.display import display, Markdown

display(Markdown(f"$$\n{latex_code}\n$$"))
$$ \begin{array}{lccc} & u & v & w \\ \hline z1 & 0.065778 & 0.050730 & 0.045735 \\ z2 & 0.041295 & 0.037394 & 0.031453 \\ z5 & 0.029445 & 0.028835 & 0.024260 \\ z10 & 0.018262 & 0.018303 & 0.014911 \\ z15 & 0.012479 & 0.012840 & 0.011146 \\ z20 & 0.008523 & 0.008437 & 0.006734 \\ z30 & 0.002335 & 0.001984 & 0.001727 \\ \end{array} $$
In [76]:
from IPython.display import display, Markdown

display(Markdown(f"$$\n{latex_code}\n$$"))
$$ \begin{array}{lccc} & u & v & w \\ \hline z1 & 0.065778 & 0.050730 & 0.045735 \\ z2 & 0.041295 & 0.037394 & 0.031453 \\ z5 & 0.029445 & 0.028835 & 0.024260 \\ z10 & 0.018262 & 0.018303 & 0.014911 \\ z15 & 0.012479 & 0.012840 & 0.011146 \\ z20 & 0.008523 & 0.008437 & 0.006734 \\ z30 & 0.002335 & 0.001984 & 0.001727 \\ \end{array} $$

The cospectrum budget model

Parameters:
  • A = 4.5
  • Γ[z] = dU/dz
  • β0 = g / θv
  • θv = mean potential temperature in Kelvin at height z; potential temperature = temperature in Kelvin * 0.0098 * z
  • Note: These are spectra as a function of wavenumber k.
  • k = 2Ï€ / λ
  • λ = wavelength = mean u wind speed / frequency
  • Ï„ = ε_u^-1/3 * k^-2/3
In [509]:
import matplotlib.image as mpimg

# Load and display the image
img = mpimg.imread("cospectrumModel.png.png")
plt.imshow(img)
plt.axis("off")  # Hide axes
plt.show()

notes:

  • i'm using the smoothed spectra
  • theta v is averaged in $\beta_0$ but not in the cospectra, since it is needed as a time series

let's take a look at the comparison data-model for $F_{uw}$ at hour 20 for each height level

In [ ]:
# evaluate the model and save some intermediate steps (e.g. the model in its simplified form)
ks, cospModel, modelNOFut, Ewws, Futs = useful.getCospectrumModel_allHeights(df_matrix_dict, hour= "h20")
In [78]:
# Fuw from data
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
Fwus = np.empty(7, dtype=object)
for i in range(7):  
    freqs, cospectrumSmooth = useful.CalcCospectrum(df_matrix_dict, hour="h20", height = heights[i], quantity1='w', 
               quantity2='u', movingAve=True, windowsize=50, highfreq = 2*10**-2)                                      
    Fwus[i] = cospectrumSmooth

in the plots we will filter out negative values for readability

In [95]:
# function to plot the results
def plotModel(altezza):
    valid_mask = ~np.isnan(ks[altezza]) & ~np.isinf(ks[altezza]) & \
                 ~np.isnan(Fwus[altezza]) & ~np.isinf(Fwus[altezza]) & \
                 ~np.isnan(cospModel[altezza]) & ~np.isinf(cospModel[altezza]) & \
                 ~np.isnan(Ewws[altezza]) & ~np.isinf(Ewws[altezza]) & \
                 ~np.isnan(Futs[altezza]) & ~np.isinf(Futs[altezza])

    # Filter y-values for the left plot (positive values only)
    left_mask_Fwus = Fwus[altezza] < 0
    left_mask_cospModel = cospModel[altezza] > 0
    left_mask_Ewws = Ewws[altezza] > 0
    left_mask_Futs = Futs[altezza] > 0
    left_mask_modelNOFut = modelNOFut[altezza] > 0

    # Set up a figure for just the left plot
    fig, ax1 = plt.subplots(figsize=(10, 4))

    # Plot the original values on the left subplot with fine lines (only positive y-values)
    ax1.loglog(ks[altezza][valid_mask & left_mask_Fwus], -ks[altezza][valid_mask & left_mask_Fwus] * Fwus[altezza][valid_mask & left_mask_Fwus], color="black", linestyle="-", linewidth=0.8, label="data")
    ax1.loglog(ks[altezza][valid_mask & left_mask_cospModel], ks[altezza][valid_mask & left_mask_cospModel] * cospModel[altezza][valid_mask & left_mask_cospModel], linestyle="-", linewidth=0.8, label="model")
    ax1.loglog(ks[altezza][valid_mask & left_mask_Ewws], ks[altezza][valid_mask & left_mask_Ewws] * Ewws[altezza][valid_mask & left_mask_Ewws], linestyle="-", linewidth=0.8, label="Eww")
    ax1.loglog(ks[altezza][valid_mask & left_mask_Futs], ks[altezza][valid_mask & left_mask_Futs] * Futs[altezza][valid_mask & left_mask_Futs], linestyle="-", linewidth=0.8, label="Fut")
    ax1.loglog(ks[altezza][valid_mask & left_mask_modelNOFut], ks[altezza][valid_mask & left_mask_modelNOFut] * modelNOFut[altezza][valid_mask & left_mask_modelNOFut], linestyle="-", linewidth=0.8, label="modelNOfut")
    
    ax1.set_xlabel(r"$k_x \, [m^{-1}]$")
    ax1.set_ylabel(r"$k_x \cdot F_{uw} \, [m/s^2]$")
    ax1.legend(fontsize=8)
    ax1.grid(True, which="both", linestyle="--", linewidth=0.5)

    # Display the plot
    plt.suptitle(f"Comparison h={heights[altezza]}")
    plt.tight_layout(rect=[0, 0, 1, 0.95])  # Adjust layout to fit title
    plt.show()
In [96]:
plotModel(0)
In [97]:
plotModel(1)
In [98]:
plotModel(2)
In [99]:
plotModel(3)
In [100]:
plotModel(4)
In [101]:
plotModel(5)
In [102]:
plotModel(6)

The behaviour of the model seems to worsen when F_uT is considered. For this reason we will study now $\tau(k)$ with the two models, complete and simplified. Be aware: the form of $\tau(k)$ we hypothesized in the model ($k^{-2/3}$ for each $k$) is a very simple one.

In [116]:
# Load and display the image
img = mpimg.imread("tauofk.png")
plt.imshow(img)
plt.axis("off")  # Hide axes
plt.show()
In [ ]:
ks, taus , tausSimple = useful.getTau_allHeights(df_matrix_dict, hour="h20")
In [108]:
# Plotting
plt.figure(figsize=(8,4))

for i in range(7):
    # Create a mask to filter out negative values in ks[i] and taus[i]
    positive_mask = (taus[i] > 0)
    
    # Plot only positive values
    plt.loglog(ks[i][positive_mask], taus[i][positive_mask], linestyle="-", linewidth=1.1, label=f"{heights[i]}")

plt.xlabel("$k_x \, [m^{-1}]$")
plt.ylabel(r"$\tau \, [s]$")
plt.title("Complete model, h20")
plt.legend()
plt.grid(True, which="both", linestyle="--", linewidth=0.5)
plt.show()
In [109]:
# Plotting
plt.figure(figsize=(8,4))

for i in range(7):
    # Create a mask to filter out negative values in ks[i] and taus[i]
    positive_mask = (tausSimple[i] > 0)
    
    # Plot only positive values
    plt.loglog(ks[i][positive_mask], tausSimple[i][positive_mask], linestyle="-", linewidth=1.1, label=f"{heights[i]}")

plt.xlabel("$k_x \, [m^{-1}]$")
plt.ylabel(r"$\tau \, [s]$")
plt.title("Simplified model, h20")
plt.legend()
plt.grid(True, which="both", linestyle="--", linewidth=0.5)
plt.show()

Extra

analysis were conducted on the hour 5 too, showing even more troubles with the $F_{uT}$ part. at the moment they are not included for brevity.

In [ ]: